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Relative motion in the Cartesian or overset framework causes certain spatial nodes 
to move in and out of the physical domain as they are dynamically blanked by moving 
solid bodies. This poses a problem for the conventional Time-Spectral approach, which 
expands the solution at every spatial node into a Fourier series spanning the period of 
motion. The proposed extension to the Time-Spectral method treats unblanked nodes 
in the conventional manner but expands the solution at dynamically blanked nodes in a 
basis of barycentric rational polynomials spanning partitions of contiguously defined tem- 
poral intervals. Rational polynomials avoid Runge’s phenomenon on the equidistant time 
samples of these sub-periodic intervals. Fourier- and rational polynomial-based differenti- 
ation operators are used in tandem to provide a consistent hybrid Time-Spectral overset 
scheme capable of handling relative motion. The hybrid scheme is tested with a linear 
model problem and implemented within NASA’s OVERFLOW Reynolds-averaged Navier- 
Stokes (RANS) solver. The hybrid Time-Spectral solver is then applied to inviscid and 
turbulent RANS cases of plunging and pitching airfoils and compared to time-accurate and 
experimental data. A limiter was applied in the turbulent case to avoid undershoots in 
the undamped turbulent eddy viscosity while maintaining accuracy. The hybrid scheme 
matches the performance of the conventional Time-Spectral method and converges to the 
time-accurate results with increased temporal resolution. 

I. Introduction 

F ORCED periodic flows arise in a broad range of aerodynamic applications such as rotorcraft, turboma- 
chinery, and flapping wing configurations. Standard practice involves solving the unsteady flow equations 
forward in time until the initial transient exits the domain and a statistically stationary flow is achieved. 
It is often necessary to simulate through several periods of motion to remove the initial transient, making 
unsteady design optimization prohibitively expensive for most realistic problems. An effort to reduce the 
computational cost of these calculations led to the development of the Harmonic Balance method [1,2] which 
capitalizes on the periodic nature of the solution. This approach exploits the fact that forced temporally 
periodic flow, while varying in the time domain, is invariant in the frequency domain. Expanding the tem- 
poral variation at each spatial node into a Fourier series transforms the unsteady governing equations into 
a coupled set of steady equations in integer harmonics that can be tackled with the acceleration techniques 
afforded to steady-state flow solvers. Other similar approaches, such as the Nonlinear Frequency Domain 
[3,4,5], Reduced Frequency [6], and Time-Spectral [7,8,9] methods, were developed shortly thereafter. Ad- 
ditionally, adjoint-based optimization techniques [10, 11] can be applied to provide the ability to perform 
design optimization without resorting to costly unsteady adjoint methods. Frequency-adaptive methods 
[12, 13, 14] can offer even greater efficiency by refining the the number of modes resolved at each grid point 
to the frequency content in its solution. 
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The Fourier temporal basis functions imply spectral convergence as the number of harmonic modes, 
and correspondingly, the number of time samples, N, increases. Some approaches solve the equations in 
the frequency domain directly. Others transform the equations back into the time domain, potentially 
simplifying the process of augmenting this capability to existing solvers. However, each approach harnesses 
the underlying steady solution in the frequency domain. These temporal projection methods will herein be 
collectively referred to as Time-Spectral methods. 

Time-Spectral methods have demonstrated marked success in reducing the computational costs associated 
with simulating periodic forced flows, but have yet to be generally applied to overset or Cartesian solvers for 
relative motion with dynamic hole-cutting. Overset and Cartesian grid methodologies are versatile techniques 
capable of handling complex geometry configurations with relative motion between components, and are 
commonly used for practical engineering applications. The combination of the Time-Spectral approach 
with this general capability potentially provides an enabling new design and analysis tool. In an arbitrary 
moving-body scenario for these approaches, a Lagrangian body-fixed mesh moves through a fixed Eulerian 
mesh. Eulerian mesh points interior to the solid body are removed ( cut or blanked ), leaving a hole in the 
Eulerian mesh. These grid points undergo dynamic blanking and do not have a complete set of time samples, 
preventing direct application of the Time-Spectral approach. Murman [6] demonstrated the Time-Spectral 
method for a Cartesian solver with rigid domain motion, wherein the hole cutting remained fixed. Similarly, 
Custer et al. [15, 16] used the NASA overset OVERFLOW solver with static hole-cutting and interpolation. 
Both of these efforts focused on applying this method to Cartesian and overset meshes with constant blanking. 
Recently, Mavriplis et al. [17] demonstrated a qualitative method for applying the Time-Spectral approach 
to an unstructured overset solver with relative motion. 

The goal of the current work is to develop a robust and consistent method for handling relative motion 
with the Time-Spectral approach within an overset or Cartesian mesh framework, while still approaching 
the spectral convergence rate of the non-overset scheme. The paper begins in Section II with a brief outline 
of the Time-Spectral approach and the issues associated with extending it to overset solvers. Section III 
details the proposed method that enables relative-body motion and provides convergence properties of the 
new differentiation operator. A linear model problem demonstrates the viability of the approach. Section IV 
outlines the implementation of the approach within OVERFLOW. The paper concludes with application of 
the Time-Spectral OVERFLOW solver to two-dimensional simulations of both plunging and pitching airfoils 
to evaluate the accuracy and convergence properties of the proposed scheme. 

II. Application of the Time-Spectral Method within an Overset Framework 

The fundamental challenge of providing Time-Spectral capabilities to an overset solver concerns the nodes 
that dynamically move in and out of the physical domain due to the relative motion between the surface 
geometry and the background Eulerian grid(s). Projection into the frequency domain has infinite support 
in time, requiring special treatment for these nodes whose solution is undefined for a subset of the temporal 
period. Figure 1 illustrates this issue for a one-dimensional oscillating piston. Fictional solutions in Fig. 
lb corresponding to nodes a, b and c in Fig. la demonstrate dynamic hole cutting by the motion of the 
piston relative to the background Eulerian grid. Node a is never blanked by the piston and therefore its 
solution is defined over the entire time history of the period, T. The piston blanks node b briefly and node 
c twice. Thus, node b has one sub-periodic time interval associated with it while node c has two - one each 
represented by the solid and hashed lines. The shaded regions serve to highlight the time over which each 
node lies outside the physical domain with an undefined solution. Special treatment is required for nodes b 
and c, while the standard Time-Spectral approach can be applied to node a, resulting in a hybrid scheme. 

For a brief synopsis of the standard Time-Spectral approach, consider a time-dependent scalar ODE 

^+R(u) = 0. (1) 

The continuous solution is expanded in the harmonic basis 

N—l 

»(<>=£ u k <t>k (t) (2) 

k=0 

from data at equispaced time samples, Uj = u(tj) = u (j At), where At = such that 

u = $ _1 u (3) 
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(a) Piston Motion (b) Solution Time Histories 

Figure 1: Figurative piston trajectory and solution histories, (a) The piston “cuts” nodes b and c as 
it oscillates over the background grid, (b) Shaded regions represent blanked intervals through which the 
solution is undefined for nodes of the corresponding color. 


with 


*jk = = e*** ( 4 ) 

serving as the transformation operator. The differentiation operator is applied in the frequency domain, i.e. 
= icukuk, and the result is projected back into the time domain leading to 

$£»$" 1 u + R(u) = 0 (5) 


where D is the diagonal Fourier differentiation operator, djk = iujkdjk- Direct assembly of the transformed 
differentiation operator, V = $D$ _1 , is derived in [7,8]. Iterating Eq. 5 in pseudotime, r, provides a 
steady-state solution procedure to an underlying unsteady, yet periodic, problem. a 

-f-u + Du + R(u) = 0 (6) 

dr 

The majority of spatial nodes in an overset or Cartesian simulation possess complete time histories and are 
thus treated with this standard procedure, but an accurate and robust treatment is still required for the 
complement of nodes that undergo dynamic blanking. 

One approach is to fill the blanked nodes with solution data via spatial averaging and proceed with the 
standard Time-Spectral method (cf. [17]). While attractive for its simplicity, this approach is inconsistent. 
Non-physical information provided by an alternative governing equation (spatial smoothing is governed by 
Laplace’s equation) is propagated into the physical domain via the infinite support of the Fourier basis. To 
demonstrate this inconsistency, see Fig. 2 which shows the solution of the one-dimensional linear advection 
equation 


du 

dt 


+ a 


du 

dx 


= 0 


( 7 ) 


given u(0,t) = — sin (uiat) for x £ [0,1]. The analytic solution is u(x,t) = sin(w(x — at)). A gap is defined 
in x to simulate blanking and the solution is linearly interpolated through the gap at every step. The 
interpolation error across the spatial gap (Fig. 2a) at one time-sample corrupts the solution throughout the 

a Alternatively, the residual operator can also be projected into the frequency domain, and the resulting equations solved 
iteratively in the frequency domain [3,6]. This approach takes advantage of the Fast Fourier Transform which is an O (ATlog N) 
operation as opposed to the O (AT 2 ) matrix operation in Eq. 5. However, for the small N typically encountered in practical 3D 
Time-Spectral applications (0(10— 100)) the difference is minimal and the time domain formulation allows for straightforward 
implementation within a steady-state solver with the addition of the temporal differentiation source term. 
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period (Fig. 2b) via the infinite support of the Fourier basis. While the resulting solution may be smooth, 
it is inaccurate and dependent upon the arbitrary averaging of the solution through the gap. Consequently, 
spatial interpolation inhibits the desired temporal convergence (Fig. 2c). If instead the node located at 
x = 0.5 were blanked, the error in the solution would be much lower because the solution at this location 
can be more accurately approximated by linear interpolation between its neighbors. 




(b) Exact and computed solutions. 


10 1 F 



log(JV) 


(c) Convergence of analytic error, e = \u — u e x |, with N 

Figure 2: Spatial smoothing Time-Spectral approach with a blanked node at x = 0.7 (a) The solution at 
the blanked node is interpolated from its neighbors, (b) The solution resulting from spatial averaging may 
be smooth but is inconsistent, as the arbitrarily filled value at the blanked node corrupts its complete time 
history and, by extension, the entire space-time domain. Severity of the corruption may be masked by 
fortuitous interpolation but it is highly dependent on the underlying solution and the gap through which it 
is filled, (c) Spatial averaging of the blanked point introduces an error that inhibits temporal convergence 
as the number of temporal modes, N , is increased. 
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III. Hybrid Fourier-Rational Time-Spectral Approach 


To develop a consistent scheme without the issues associated with spatial averaging presented in the 
previous section, the current approach partitions the time histories for spatial nodes exhibiting dynamic 
blanking into intervals of consecutively unblanked time samples. Global temporal support is thus abandoned, 
however, this approach is consistent with the physics of disjoint domains separated by an impermeable 
boundary. It is important to note that these temporal partitions must still be evaluated on a uniform lattice 
to avoid manipulating the spatial operator. Thus the problem is now reduced to finding a robust and efficient 
method for evaluating the temporal derivatives of non-periodic functions on a uniform lattice. 

Several approaches to achieve this goal were evaluated. One option is to compute derivatives with 
finite differences or a localized differentiable basis such as wavelets. However, compact bases offer limited 
accuracy [18] and wavelet differentiation requires special treatment at non-periodic interval boundaries [19]. 
Another option is temporal extrapolation that populates the undefined portion of the period using data 
from the physical portion of the time signal. Examples of this approach, such as Fourier continuation, which 
extrapolates a non-periodic function into a periodic function on a larger domain [20,21,22], and compressed 
sampling, which requires L\ minimization to solve an underdetermined system [23], are both too costly and 
lack the requisite robustness for this application. 

The current approach projects the solution of the unblanked nodes within a single partition onto a 
global basis spanning the same interval. The constraint of evenly-spaced temporal collocation points elimi- 
nates traditional optimally-accurate orthogonal polynomials ( e.g . Chebyshev) due to Runge’s phenomenon 
at the endpoints. A least-squares projection of orthogonal polynomials is more stable, but results in an 
overdetermined system whose projection does not interpolate the solution data. Barycentric rational poly- 
nomials provide a viable alternative [24] and Bos et al. [26] showed that they offer superior interpolation 
and differentiation properties on equidistant nodes over conventional orthogonal polynomials. A rational 
polynomial-based differentiation operator will therefore be used in this hybrid Time-Spectral approach. An 
alternate basis may prove more successful in the future so the current approach will be kept general enough 
to use any projection operator, contingent on the availability of its analytic differentiation operator. In the 
proposed scheme, solutions across each non-periodic interval are projected onto a basis of rational polynomi- 
als and differentiated accordingly. Fully periodic intervals are still projected and differentiated in the Fourier 
basis, resulting in a hybrid approach that employs the optimal basis available. 

Floater and Hormann [25] introduced a set of barycentric rational polynomials with high rates of ap- 
proximation convergence and no poles. Additional efforts to explore rational polynomials and their utility 
as a pseudospectral basis for spectral collocation methods include but are not limited to [27,28,29,30]. 
Baltensperger et al. [31] define the rational interpolant, r(x ), in barycentric form 

N 

E 

r(x) = ^ ( 8 ) 

E Wk 
X-Xk 

k—0 


and its corresponding differentiation operator, D, such that Ar = Dr 


Djk — 


Wk 1 
Wj ( Xj—Xk ) 

N 

- E D ri 

i=0,i^k 


if j ^ k 
if j = k 


where rk = r(xk)- The interpolation can be reformulated as a weighted sum of nodal basis functions, 4>{x), 
with coefficients, /, equal to the function value at each node, fk = /( Xk)- 


N 

f ( x ) = Y fk4>k(x), with <t> k (x) 

k—0 


Wk 


X — Xk 


N 

ip Wj 

E— / x—Xj 
3=0 3 


(9) 


Barycentric rational polynomials nodally interpolate the solution data resulting in an identity transformation 
operator, = (f>k(xj ) = 6jk, implying discrete orthogonality. It follows that the transformed differentiation 
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operator is identical to D: V = l*” 1 = D. While similar in form to the Lagrange interpolant, a 

key distinction of the barycentric rational polynomials is in how the weights, Wk, are defined. Lagrange 
polynomials are constrained to pass through TV + 1 points as a polynomial of order TV. In contrast, while 
the rational polynomials interpolate the data, they are not constrained to do so as a polynomial of order 
TV. This relaxation aids in preventing spurious oscillations at the endpoints of equispaced nodes associated 
with the traditional Lagrange or Chebyshev polynomial bases while retaining powerful interpolation and 
differentiation properties. The weights, Wk, control the accuracy and stability characteristics of the rational 
interpolant. Floater and Hormann [25] derived weights that provide an approximation order of d while 
avoiding poles. Weights guaranteeing an absence of poles 


w k 


i-\-d 

n 

ieJk j=i,j^k 


1 

\x k - Xj I 


( 10 ) 


are defined for TV + 1 samples where 


It ’■= {0, 1, ..., TV — d} and J a := {* £ It '■ ot — d < i < a} (11) 

for the desired approximation order d. Readers are directed to the derivation in [25] for additional details. 
Figure 3 illustrates three basis functions corresponding to the first, second and fourth nodes with approxi- 
mation order, d = 3, over 8 sample points (TV = 7) with the values highlighted at the nodal points. Note the 
interpolation property of the basis, 4>k(xj) = Sjk- 



x 

Figure 3: Barycentric rational polynomial basis functions <j>o{x), <f> i(x) and <p r i{x) corresponding to nodes 0, 
1 and 3, respectively, for approximation order, d = 3, and 8 equispaced sample points. 

Differentiation performance of the rational polynomial basis is compared to that of third-order finite- 
differences, the Fourier differential operator (generally optimal for periodic functions), differentiation using 
least-squares orthogonal polynomials on evenly-spaced nodes (modes M < 2y/N) and the Chebyshev differ- 
ential operator (generally optimal for non-periodic functions) on Chebyshev nodes. An even-odd harmonic 
function and Runge’s function (designed to expose numerical instabilities) are differentiated by each of the 
aforementioned methods. Convergence of differentiation error, e = ||^ — X>/|| 2 , with TV is plotted in Fig. 
4 for the harmonic function and Fig. 5 for Runge’s function. The approximation order for the rationals is 
selected as d = min ([^J ,d max ) with d max = 10 for this case, but optimal selection is problem dependent 
and an area of continuing research [29,32]. The rational polynomials approach spectral convergence for both 
the periodic and non-periodic problems. Fourier differentiation exhibits roundoff error for the harmonic 
function with increased TV and offers poor convergence for Runge’s function. 
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log(JV) 



(a) Analytic function and derivative 


(b) Convergence of differentiation operators 


Figure 4: Convergence of differentiated harmonic function, /( x) = cos(wa;) + sin(w:r) 



x 


1 Third-order 
■ Spectral 
' Rational 
1 Least Squares 
1 Chebyshev 
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log (iV) 


(a) Analytic function and derivative 


(b) Convergence of differentiation operators 


Figure 5: Convergence of differentiated Runge’s function, f(x) = 


l 

l+25a; 2 


To demonstrate the viability of the proposed approach, the same one-dimensional linear advection case 
(Eq. 7) shown in Sec. II is used to test the hybrid Time-Spectral strategy. Instead of the spatial interpolation 
used previously, the barycentric rational polynomial temporal differentiation operator, T > r , is used on the 
spatial nodes blanked over some duration of the period and the Fourier differentiation operator, is used 
on nodes with a periodic interval. Information propagates from left to right with wave speed a = 1 and thus 
analytic boundary conditions are supplied both at * = 0 and along the right edge of the blanked region. The 
domain is discretized with N x = 101 points in space and N = 21 time samples. The solution is initialized 
uniformly to zero outside the boundary nodes, and third-order finite-differences are used to approximate 
spatial derivatives to minimize errors associated with the spatial discretization. The equations are solved 
implicitly and converged to machine precision. The space-time solution, visualized in Fig. 6, is compared to 
the analytic solution and the solution of the unblanked case. 
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(a) Space-time solution 




X 

(b) Solution and blanked interval at fixed time t = 0.4286 

Figure 6: Solution for linear advection test problem using the hybrid Fourier-rational polynomial Time- 
Spectral approach. The rational polynomial differentiation operator, T> r , is applied to nodes in the central 
region of the domain that undergo blanking and the Fourier differentiation operator, T> f, is applied elsewhere. 

The maximum error between the solution computed with the hybrid scheme and the analytic solution is 
7.757 x 10~ 3 . This value is just slightly larger than the maximum error between the unblanked case using 
only the Fourier differentiation operator and the analytic solution, 7.717 x 10 -3 . The largest difference 
between the two solutions is 5.485 x 10 -3 suggesting that the Fourier-rational hybrid approach essentially 
produces the same result as the conventional Tinre-Spectral method in the active region. 
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IV. Implementation of the Hybrid Time-Spectral Approach in OVERFLOW 


The Time-Spectral approach requires no modification to the spatial operator, making the scheme straight- 
forward to implement and maintain within an existing code. The current work leverages the OVERFLOW 
overset flow solver [33], which provides a mature code base for solving the three-dimensional compressible 
Reynolds-averaged Navier-Stokes (RANS) equations using a structured finite-difference scheme. OVER- 
FLOW employs an approximate-factorization algorithm for the implicit solution of the spatial derivatives, 
and the Time-Spectral operator is easily integrated within this framework by the addition of a similar factored 
operator for the time derivative. This approximate factorization avoids forming and solving a N x x Nq x N 
system of equations coupled across time and space. Writing the full Time-Spectral approximate-factorization 
difference equation in “delta” form gives 


[I + A tA x ] [I + AtA v \ [I + A tA z \ [I + AtV } A Q 


—Ar 


R (Q n ) + 5 (Q n ) 


( 12 ) 


where R(Q n ) is the standard nonlinear spatial residual operator, and S{Q n ) is the Time-Spectral differential 
operator applied to the state vector at time level n. A x . A y and A z are the linearized spatial operators 
corresponding to the x-, y- and z-directions. The Time-Spectral approach as applied to OVERFLOW avoids 
modifying the existing spatial residual and implicit operators because they are applied to each time sample 
sequentially within each iteration. The required modifications are limited to an additional linear solve of 
dimension N for the Time-Spectral factored operator, and computation of the source term at every spatial 
grid point. In other words, time is treated as an additional spatial dimension when solving for the steady-state 
solution in the combined space-time domain. 

The computational memory requirements for a Time-Spectral solution are greater than its time-accurate 
analog because the solution and source term must be stored for each time sample. However, careful im- 
plementation limits the storage requirements and allows for a significant number of temporal modes even 
in complex three-dimensional problems. Ignoring the trivial storage and processing associated with hole 
cutting, the current implementation in OVERFLOW stores a copy of the state vector Q n and Time-Spectral 
source term S(Q n ) at each time level, along with an additional working copy of the grid and metrics. Thus, 
the total additional storage is 2 N x NqN + 15 N%, where N x and N represent the number of spatial and 
temporal nodes, respectively, and Nq is the number of conserved solution variables. An additional cost 
of the approximate-factored matrix for the temporal derivative is also accounted for, depending upon the 
relative size of N x and N. Using a target problem size of N x = 125M grid points, which is representative 
of a three-dimensional rotorcraft or turbomachinery simulation, and 128 8-core Intel Nehalem nodes of the 
NASA Advanced Supercomputing Pleiades supercomputer, Fig. 7 shows the variation in total memory with 
number of temporal inodes for our target problem. As seen, even limiting to 128 compute nodes, with 24 
GB of memory per node, there is still sufficient memory for up to 200 temporal modes, and larger problems 
similarly benefit from the low memory requirements. 



Figure 7: Memory usage estimate (GB) versus N for 128 Intel Nehalem Nodes on the NASA Advanced 
Supercomputing (NAS) Pleiades supercomputer for a sample problem of 125 million grid points. 
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V. Numerical Results 


Two-dimensional computational studies were performed to demonstrate the ability of the hybrid Time- 
Spectral scheme to handle dynamic hole-cutting associated with overset relative motion. First, a large- 
amplitude inviscid plunging airfoil example verifies its performance in simulations involving large regions of 
dynamically blanked points. Next, an inviscid oscillating airfoil case corresponding to the AGARD 702 test 
is simulated. Finally, the turbulent case of the AGARD 702 oscillating airfoil is considered to examine how 
the added complexity of turbulence affects the hybrid scheme and the Time-Spectral approach in general. 

Results for both rigid and relative motion were computed by the OVERFLOW Time-Spectral solver 
providing a direct comparison between the time-accurate approach, the conventional Time-Spectral approach, 
and the proposed method. These problems serve as computationally tractable demonstrations of capability 
and are representative of more complex three-dimensional applications. 

Standard OVERFLOW inputs of second-order central differencing for the convective flux terms, first- and 
third-order artificial dissipation and the ARC3D Beam- Warming tridiagonal implicit spatial operator were 
employed. The second-order backwards difference operator (BDF2) is used to advance the time-accurate 
simulations. The rational polynomial order, d, is defined for all OVERFLOW Time-Spectral simulations as 
d = min ([^w^j , d max ) , with d ma x = 6 and N collocation points (not N + 1). 

Inviscid Plunging NACA 0012 Airfoil 

A large-amplitude inviscid plunging airfoil test case provides a meaningful demonstration of the hybrid Time- 
Spectral scheme as it results in a large number of dynamically blanked nodes. Plunging amplitude is chosen 
as half the chord length resulting in the following definition of plunging motion, h ( t ) , 

h (t) = ^ sin (kt) (13) 

with reduced frequency, k = 0.1627 radians per non-dimensional time unit. The free stream Mach number 
is chosen as = 0.5 to maintain subsonic flow throughout the domain. 


I 
I 

(a) Differentiation Operator Locations (b) Length of Temporal Interval at t = 0 

Figure 8: Plunging NACA 0012 airfoil, (a) The Fourier differentiation operator is applied in the yellow region 
where spatial nodes have complete time histories. The barycentric rational polynomial-based differentiation 
operator is applied to dynamically blanked spatial points shown in pink. Labels a and b locate nodes whose 
solutions are plotted in Figs. 9 and 10, respectively, (b) Length of time interval associated with each spatial 
node at the t = 0 temporal collocation point for the N = 9 case. The white region contains nodes blanked 
by the airfoil grid at the current time level. 

Figure 8a shows the portion of the domain that undergoes dynamic blanking for the plunging airfoil for 
the case of N = 9 time samples. The airfoil sweeps through the pink-colored region resulting in blanked 
spatial nodes for some subset of the period of motion. The temporal collocation points of these spatial nodes 
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are therefore partitioned into intervals spanning fewer than TV = 9 time samples and differentiated with the 
rational polynomial operator of the appropriate rank. Nodes in the yellow-colored region remain unblanked 
for the duration of the period and are therefore differentiated by the Fourier operator. The length of the time 
interval associated with each dynamically blanked point varies and is plotted in Fig. 8b for the temporal 
collocation point at t = 0. The white region shows the nodes blanked by the airfoil grid in its neutral position 
at t = 0. Nodes colored red have complete time histories and correspond to the yellow-colored nodes in Fig. 
8a. The solution at these nodes can be expanded in a Fourier series and are differentiated by the Fourier 
operator. Nodes colored yellow, green and blue are blanked for a portion of the period and only have defined 
solutions at a subset of the time samples. For example, nodes colored green are included in temporal intervals 
spanning 7 temporal collocation points including the t = 0 time sample. The distribution of interval length 
changes as the number of modes is increased but this example is representative of the general case. 

Figures 9 and 10 track the time history of the conserved solution variables at nodes a and 6, respectively, 
highlighted in Figs. 8a. Node a is located near the boundary between the blanked and unblanked regions 
and is therefore strongly influenced by the Fourier basis. Node a is blanked briefly as the airfoil sweeps 
upwards and again as it moves downwards resulting in two contiguous time intervals in the time-accurate 
case. Figure 9 plots the time-accurate and Time-Spectral solutions at node a. In the case of TV = 5 Time- 
Spectral modes, node a maintains a complete time history and is therefore differentiated by the Fourier-based 
operator. All but one temporal collocation point is active in the TV = 9 case resulting in a single partition 
that is differentiated by the rational polynomial-based operator spanning eight nodes. Two independent 
temporal partitions are generated for the TV = 17 case resulting in one interval spanning just two temporal 
collocation points and another spanning an additional twelve collocation points. Excellent agreement for all 
the modal cases is observed with only small deviations from the time-accurate solution. Node b located in 
the center of the blanked region is farther from Fourier-based nodes. The continuous solution at node b is 
partitioned into two intervals of similar length. Only the first temporal collocation point is blanked for the 
TV = 5 and TV = 9 Time-Spectral cases whereas the time history of the TV = 17 case is partitioned into two 
intervals, each spanning seven time samples. The time-accurate and Time-Spectral solutions at node b are 
plotted in Fig. 10. Larger discrepancies between the time-accurate and Time-Spectral solution of streamwise 
momentum, pu , are manifested for the TV = 5 and TV = 9 cases (Fig. 10b) but the solution does converge 
with increased temporal resolution evidenced by the TV = 17 case. 
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Time Accurate 



t/T 

(a) Density, p (t) 



t/T 

(b) x-Momentum, pu (t) 



t/T 

(c) z-Momentum, pw (t) 

Figure 9: Time-accurate and Time-Spectral solutions of p , pu and pw for the plunging NACA 0012 airfoil 
at node a in Fig. 8a. Values from the Time-Spectral solutions using N G {5, 9, 17} time samples are plotted 
at their respective collocation points against the time-accurate solutions. 
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(b) ^-Momentum, pu ( t ) 



t/T 

(c) z-Momentum, pw (t) 

Figure 10: Time-accurate and Time-Spectral solutions of p , pu and pw for the plunging NACA 0012 airfoil 
at node b in Fig. 8a. Values from the Time-Spectral solutions using N £ {5,9, 17} time samples are plotted 
at their respective collocation points against the time-accurate solutions. 
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(a) N = 5 (b) N = 9 (c) N = 17 



■ 

Legend 


Figure 11: Subsonic plunging NACA 0012 airfoil. Normalized error in pressure, err p = ^ Pt ^ t Pts ^ , of the 
Time-Spectral versus time-accurate solutions at t = 0 for TV S {5,9,17}. Error is plotted on a log scale 
ranging from -2 to -7 suggesting that even the N = 5 case has less than 1% error in pressure over then entire 
flow-field at the t = 0 instance. 



h/c 

Figure 12: Subsonic plunging NACA 0012 airfoil. Convergence of relative-motion Time-Spectral projection 
of drag coefficients towards time-accurate solution for N € {3,5,9,17}. The Time-Spectral drag signals 
for the N = 9 and N = 17 cases are observed to nearly match the time-accurate solution demonstrating 
temporal convergence. 


To examine solution quality on a larger scale the global error in pressure is computed taking the relative- 
motion time-accurate pressure, pta, as exact. Normalized error in pressure, err p , is visualized in Fig. 11 for 
N £ {5, 9, 17} at time t = 0 which is the only common temporal collocation point among the three examples. 
The maximum normalized error in pressure occurs in the N = 5 flow held but is limited to roughly one part 
in 1000. The error is reduced monotonically by increasing the number of modes. This is evident in the 
drag coefficient in Fig. 12. The N = 3 and N = 5 Time-Spectral cases demonstrate differences from the 
time- accurate result, but the N = 9 case shows improvement and the A = 17 case matches the time- accurate 
solution nearly identically. 

Transonic AGARD 702 Pitching Airfoil 

The Time-Spectral formulation is next compared to time-accurate simulations and experimental data corre- 
sponding to the AGARD 702 3E3 oscillatory pitch test case CT5 [34]. Both inviscid and viscous turbulent 
simulations are performed on this transonic, M, ^ = 0.755, case. A pitching airfoil allows for simulations 
with both rigid- and relative-body motion so that the temporal discretization can be isolated and compared 
directly. The airfoil pitches about its quarter chord with incidence, a (t) , defined as 

a (t) = 0.016° + 2.51° sin (kt) (14) 
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with reduced frequency k = 0.1627 radians per non-dimensional time unit and free stream Mach number 
Moo = 0.755. A hierarchy of Time-Spectral simulations are computed using an increasing number of modes 
(N € {3,5,9,17,33}) to investigate the performance of the scheme with increased temporal resolution and 
ensure that spurious higher-frequency modes are not generated and propagated into the solution. The near- 
body grid and its overlap region are chosen such that the transient shock spans the grid interface. Figure 
13 shows the supersonic region and the overset-grid system at two of the temporal collocation points for the 
N = 5 case. 



Figure 13: AGARD 702 inviscid pitching airfoil at = 0.755. The supersonic region is plotted on the 
grid-system in the vicinity of the airfoil for two time samples of the N = 5 case demonstrating the shock 
spanning the near-body to off-body grid interface. Subsonic regions shown in black. Every second grid point 
has been removed in the figure for clarity. 


Inviscid Case 

A 241 x 30 O-mesh near-body grid is embedded in a 341 x 261 rectilinear background grid that stretches 
approximately 100 chord lengths to the far field in both directions. In the case of rigid-body motion both 
the near-body and background grids rotate together. This provides constant hole-cutting permitting the 
use of the conventional Time-Spectral approach using the Fourier-based differentiation operator. In the 
relative-motion case, the near-body grid rotates independently of the background grid introducing dynamic 
hole-cutting. The unsteady residual was converged to machine precision at each physical time step in the 
time-accurate case. The complete space-time residual for the Time-Spectral cases also converges to machine 
precision (Fig. 14). The convergence rate is independent of the number of Time-Spectral modes in this 
case. 



Iterations x 10 4 


Figure 14: AGARD 702 inviscid pitching airfoil. Iterative convergence of space-time residual for relative- 
body motion for N £ {3, 5, 9, 17, 33}. The sum of both the near- and off-body grid residuals normalized by 
the total number of points in the domain is plotted versus iteration number. 
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Pitching moment and drag coefficient polars are plotted in Figs. 15 and 16, respectively. The rigid- and 
relative-motion Time-Spectral solutions and the time-accurate case, all converge to the identical solution. 
Reasonable agreement with the experimental data is achieved but most importantly, significant temporal 
resolution convergence of the hybrid Time-Spectral solutions is achieved as the number of sample points 
is increased from 3 to 33 allowing the solution to resolve increasingly higher-frequency content inherent in 
the unsteady solution and its output functionals. Note that 33 modes are required to resolve the kink in 
the pitching moment and drag signals at a ± 2° . Root-mean-square error between the unsteady force and 
pitching moment coefficients of the rigid- and relative-motion case isolates the quantitative effect of the 
hybrid differentiation operator on the Time-Spectral solution and confirms the convergence (Table 1). 


N 


c d 

C m 

3 

1.96e-03 

3.64e-04 

6.24e-04 

5 

3.50e-04 

3.83e-05 

2.01e-04 

9 

7.83e-04 

1.24e-05 

4.22e-05 

17 

6.98e-05 

1.03e-05 

1.68e-05 

33 

4.50e-05 

4.08e-06 

1.39e-05 


Table 1: Root-mean-square error in force and pitching moment coefficients between projections of rigid- and 
relative-motion Time-Spectral solutions. 


The frequency components, amplitude and phase, of the time-accurate pitching moment coefficient are 
plotted in Fig. 17 against those computed by the hybrid Time-Spectral scheme. Figure 17a reveals that 
the N = 3 case fails to resolve even the mean amplitude and the N = 5 case poorly resolves the either the 
amplitude or phase corresponding to the second harmonic. Increasing the number of solution modes permits 
resolution of higher frequency content, and improves resolution of the lower frequencies. Importantly, no 
artificial high-frequency amplitude or phase content corrupts the solution as a result of the hybrid Time- 
Spectral scheme so that increasing the number of solution modes results in solutions that mirror the frequency 
content of the time-accurate simulation. 
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(a) N = 3 




(c) N = 33 

Figure 15: AGARD 702 inviscid pitching airfoil. Time Spectral versus Time Accurate pitching moment 
coefficients for N £ {3, 9, 33}. Ten periods of the time-accurate solution are plotted in red from steady-state 
startup. Blue squares locate the pitching moment coefficient values at the Time-Spectral collocation points 
for relative-body motion. Relative-body pitching moment coefficients computed from an interpolation of the 
Time-Spectral solution to 201 points shown with the blue-hashed line. Green diamonds locate the pitching 
moment coefficient values at the Time-Spectral collocation points for rigid-body motion. Rigid-body pitching 
moment coefficients computed from an interpolation of the Time-Spectral solution to 201 points shown with 
the green-hashed line. Experimental data from the AGARD 702 Report is plotted with black dots. 
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a 

(a) N = 3 




a 

(c) N = 33 

Figure 16: AGARD 702 inviscid pitching airfoil. Time Spectral versus Time Accurate drag coefficients for 
N £ {3,9,33}. Ten periods of the time-accurate solution are plotted in red from steady-state startup. Blue 
squares locate the drag coefficient values at the Time-Spectral collocation points for relative-body motion. 
Relative-body drag coefficients computed from an interpolation of the Time-Spectral solution to 201 points 
shown with the blue-hashed line. Green diamonds locate the drag coefficient values at the Time-Spectral 
collocation points for rigid-body motion. Rigid-body drag coefficients computed from an interpolation of the 
Time-Spectral solution to 201 points shown with the green-hashed line. Corrected experimental drag data 
from the AGARD 702 Report is not available. 
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(a) Amplitude, 5R|C m | 



(b) Phase, Qf|Crrc| 

Figure 17: AGARD 702 inviscid pitching airfoil. Spectrum of time-accurate and Time-Spectral pitching 
moment coefficients, C m = for first 20 harmonics of the relative motion case. Increasing modal 

resolution results in better agreement between output functional amplitude and phase without introducing 
spurious high-frequency content. 


Turbulent RANS Case 

Practical applications of the Time-Spectral method including rotorcraft and turbomachinery involve turbu- 
lent flows, and therefore the previous oscillatory airfoil is extended to use the Spalart-Allmaras one-equation 
turbulence model. OVERFLOW employs a loosely coupled turbulence scheme whereby the turbulent vari- 
ables are updated initially and held constant for the flow equation iteration. The code structure prevents a 
direct Time-Spectral implicit treatment for the turbulent variable without a significant overhaul. Instead, a 
semi- implicit treatment retroactively applies the implicit operator to the explicit turbulent update, A£' ra_t "5. 

[I + AtA x ] [I + At A y ] [I + AtA 2 ] Au n+ i = -At [R (P n ) + S (£>")] (15) 

( 16 ) 
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v n+ 2 is then held fixed while iterating the flow solution, Q. 

Q n+1 = Q n + AQ (v n+ i'j (17) 

The implicit update Av n+1 is computed and retroactively applied to v n to advance the solution prior to the 
subsequent time step. 

[I + AtV] AD n+1 = Au n+ 3 (18) 

(19) 

Figure 18 shows contours of the (undamped) turbulent eddy viscosity, u, in the wake region for the five 
temporal collocation points of the pitching airfoil. This illustrates the potentially large changes in v as a 
function of time in the wake region. 



Figure 18: Visualization of undamped turbulent eddy viscosity on the near-body grid and in the wake region 
of the AGARD 702 pitching airfoil at the temporal collocation points for a representative case of N = 5. 
v can vary several orders of magnitude in a short time-frame and must be limited in certain circumstances 
to avoid spurious undershoots resulting from Gibbs’ phenomenon. The figure shows every second grid point 
removed in the normal direction for clarity. 


Figure 19 plots the time history of turbulent eddy viscosity and its temporal derivative at a point 
downstream of the airfoil through which the turbulent wake passes during a portion of the oscillation. Note 
in Fig. 19a, the turbulent eddy viscosity transitions from a constant value of order zero to a value of order 
1000 over a small fraction of the period. While the function varies smoothly in the context of the small 
physical time steps associated with time-accurate iterations, it resembles a discontinuous step function in 
the frame of the much coarser resolution of Time-Spectral sample points. This leads to spurious oscillations 
in the Fourier expansion (Gibbs’ phenomenon). Oscillations in the gradient are manifested at the discrete 
sample points which can be observed in Fig. 19b resulting in large overshoots. This issue is compounded 
by the fact that undamped eddy viscosity remains close to zero outside the wake, so that modest overshoots 
in the derivative violate the positivity constraint on u. The current example uses eddy viscosity in the 
Spalart-Allmaras model, however, similar issues are encountered for turbulent kinetic energy or dissipation 
in other models. The same phenomenon could also occur in the fluid equations for high-speed cases. Note 
that this issue applies to Time-Spectral treatment of turbulence in general as it occurs downstream from the 
body on nodes with complete time histories. 

Applying a limiter to the temporal differentiation operator maintains positivity on £ without sacrificing 
accuracy. The current limiter sets the Time-Spectral source term to zero if it would drive v negative. 
Figure 19a demonstrates good agreement between the Time-Spectral and time-accurate solutions for i> at 
the collocation points. Slight ringing does occur in the constant region but v remains positive. 
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(b) Time derivative, (t) 

Figure 19: AGARD 702 turbulent pitching airfoil, (a) Turbulent solution variable, (t), and (b) its derivative, 
are plotted for a point in the wake over the period of oscillation. BDF2 and Fourier operators differenti- 
ated the continuous and discrete data, respectively. Note, Gibbs’ phenomenon in the Fourier representation 
of the signal are present in the derivative where large oscillations lead to an inaccurate representation of ^ 
(t*=t-T/ 4). 


Only the relative-motion Time-Spectral scheme is applied to the viscous case as the hybrid scheme 
already demonstrated its ability to match the results computed with the conventional Time-Spectral scheme 
in inviscid rigid-body case. The drag and pitching moment coefficients computed from the Time-Spectral 
simulations converge to those computed in unsteady mode with increasing numbers of temporal collocation 
points (Figs. 20 and 21). Note that TV = 33 modes are required to resolve the kink in the drag and pitching 
moment coefficient signals at a w ±2° (Figs. 20c and 21c), as with the inviscid case. Less agreement is 
observed with the experimental data in the pitching moment coefficient polars (Fig. 20) when compared 
to the inviscid case from the previous section but this discrepancy also occurs in the time-accurate case 
indicating suggesting a modeling issue. While experimental data is included for reference, the crucial point 
remains that the Time-Spectral results matched those computed with the time-accurate solver. 
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(a) N = 3 



(b) N = 9 



(c) N = 33 

Figure 20: AGARD 702 turbulent pitching airfoil. Time Spectral versus Time Accurate pitching moment 
coefficients for N £ {3, 9, 33}. Ten periods of the time-accurate solution are plotted in red from steady-state 
startup. Blue squares locate the pitching moment coefficient values at the Time-Spectral collocation points 
for relative-body motion. Relative-body pitching moment coefficients computed from an interpolation of the 
Time-Spectral solution to 201 points shown with the blue-hashed line. Experimental data from the AGARD 
702 Report is plotted with black dots. 
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a 

(a) N = 3 



(b) N = 9 



(c) N = 33 

Figure 21: AGARD 702 turbulent pitching airfoil. Time Spectral versus Time Accurate drag coefficients for 
N £ {3,9,33}. Ten periods of the time-accurate solution are plotted in red from steady-state startup. Blue 
squares locate the drag coefficient values at the Time-Spectral collocation points for relative-body motion. 
Relative-body drag coefficients computed from an interpolation of the Time-Spectral solution to 201 points 
shown with the blue-hashed line. Corrected experimental drag data from the AGARD 702 Report is not 
available. 
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VI. Summary 


A new hybrid Time-Spectral method has been introduced for dynamic overset and Cartesian applica- 
tions. The time samples of spatial nodes that undergo dynamic blanking are segmented into intervals of 
consecutively unblanked temporal collocation points and differentiated by an operator derived from barycen- 
tric rational polynomials. Barycentric rational polynomials offer superior interpolation and differentiation 
properties on an equispaced aperiodic mesh in contrast to traditional orthogonal polynomials. The majority 
of nodes are unblanked and treated by the standard Time-Spectral formulation. 

The hybrid scheme was outlined and evaluated with a linear model problem and the implementation of the 
hybrid Time-Spectral scheme in the NASA OVERFLOW solver was described. Inviscid and turbulent test 
cases of plunging and pitching airfoils solved using both rigid and relative motion agreed with time-accurate 
solutions. The hybrid Time-Spectral scheme performed well in the inviscid subsonic plunging airfoil case 
where its conserved solution variables and drag coefficient values converged to the time-accurate data. The 
AGARD 702 inviscid pitching airfoil case demonstrated temporal convergence in its pitching moment and 
drag signals requiring N = 33 modes to resolve the finest details of the solution. Spectrum analysis confirmed 
that no artificial high-frequency content was generated. A limiter was employed in the turbulent case to 
avoid large oscillations in the gradient of the turbulent eddy viscosity in the wake. Temporal convergence 
was similarly observed for the viscous case. The hybrid Time-Spectral scheme matched the results computed 
with both the conventional Time-Spectral approach and the time-accurate solver, establishing it as a viable 
tool for the simulation of periodically forced flows. 

Future work includes acceleration of the code through multigrid and other techniques to enable application 
to more realistic problems in rotorcraft and turbomachinery. Dynamic near-body to near-body hole-cutting 
occurs in these problems and will therefore be investigated. Optimal selection of the rational polynomial 
approximation order, d , and its maximum value, d max , remains part of future work in addition to the stability 
analysis for the barycentric rational polynomial-based differentiation operator. 
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